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Abstract 

A sequential filtering algorithm is presented for attitude and attitude-rate estima- 
tion from Global Positioning System (GPS) differential carrier phase measurements. 

A third-order, minimal-parameter method for solving the attitude matrix kinematic 
equation is used to parameterize the filter’s state, which renders the resulting estima- 
tor computationally efficient. Borrowing from tracking theory concepts, the angular 
acceleration is modeled as an exponentially autocorrelated stochastic process, thus 
avoiding the use of the uncertain spacecraft dynamic model. The new formulation 
facilitates the use of aiding vector observations in a unified filtering algorithm, which 
can enhance the method’s robustness and accuracy. Numerical examples are used to 
demonstrate the performance of the method. 

I Introduction 

A TTITUDE determination methods using Global Positioning System (GPS) signals 

have been intensively investigated in recent years. In general, these methods can 

be classified into two main classes. Point estimation algorithms (also called “deterministic” 
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algorithms), in which the GPS measurements at each time point are utilized to obtain an at- 
titude solution independently of the solutions at other time points, were introduced, among 
others, in Refs. 1, 2 and 3. Stochastic filtering algorithms, which process the measurements 
sequentially and retain the information content of past measurements, can produce better 
attitude solutions by more effectively filtering the noisy measurements. Such algorithms 
were recently introduced in Refs. 4 and 5, both of which utilized extended Kalman filtering 
to sequentially estimate the attitude from GPS carrier phase difference measurements. Both 
attitude and attitude-rate were estimated, and the filters used the nonlinear Euler equa- 
tions of motion for attitude propagation. While avoiding the traditional usage of the costly 
and unreliable gyro package, this approach rendered the resulting filters computationally 
burdensome and sensitive to inevitable modeling errors. 6 In Ref. 4 an attempt was made to 
robustify the dynamics-based filter by estimating the unknown disturbance torques, modeled 
as unknown constants. 

Although GPS-based attitude estimation methods should enjoy, in principle, the low price 
and low power consumption of state-of-the-art GPS receivers, and the general availability and 
robustness of the global positioning system, these methods are very sensitive to multipath 
effects and to the geometry of the antennae baseline configuration, and they inherently rely 
on precise knowledge of the antennae baselines in the spacecraft body frame. On the other 
hand, methods based on vector observations have reached maturity and popularity in the 
last three decades. However, as is well known, they too suffer from disadvantages, that can 
be attributed to the particular attitude sensors on wh ch they are based. Thus, while their 
readings are relatively noiseless, Sun sensors are very sensitive to Earth radiation effects, 
and axe rendered completely useless during Eclipse. Star trackers can provide accuracy on 
the order of a few arc-seconds, but are usually extremely expensive. Magnetometers always 
provide measurements of the Earth magnetic field in spacecraft flying in low Earth orbits, but 
they are sensitive to unmodeled residual magnetic fields in the spacecraft and to magnetic 
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field model imperfections and variations. 

The method presented herein is a sequential estimator for both the spacecraft attitude 
matrix and attitude-rate, which mainly uses differential GPS carrier phase measurements, 
but can also process aiding vector observations (such as low accuracy coarse Sun sensor 
measurements, or magnetic field measurements). Conceptually similar to the principle of 
complementary filtering, 7 the idea underlying this estimator is that, due to the different 
nature of these signals, the combination of both in a unified data processing algorithm can 
benefit from the relative advantages of both sensor systems, while alleviating the disadvan- 
tages of both. 

The new estimator is based on a third-order minimal-parameter method for solving the 
attitude matrix evolution equation using integrated-rate parameters (IRP). 8 Similarly to 
Refs. 5 and 4, the new estimator is a sequential filtering algorithm and not a determin- 
istic (point estimation) algorithm. However, the new algorithm differs from other works 
addressing the same problem in two main respects. First, the estimator’s propagation model 
does not utilize the nonlinear Euler equations. Instead, employing an approach borrowed 
from linear tracking theory, 9 the uncertain dynamic model of the spacecraft is abandoned, 
and the angular acceleration is modeled as a zero-mean stochastic process with exponential 
autocorrelation. [A similar, but simpler, approach was employed in the Applied Technol- 
ogy Satellite 6 (ATS-6) 10 ]. Combined with the extremely simple evolution equation of the 
integrated-rate parameters, this results in a simple, linear propagation model. Second, in 
contrast with other methods relying mainly on the attitude quaternion, the algorithm pre- 
sented herein directly estimates the attitude matrix, a natural, nonsingular attitude repre- 
sentation. Building upon the minimal, third-order integrated-rate parametrization, the new 
estimator assigns just three state variables for the parametrization of the nine-parameter 
attitude matrix, which is at the heart of its computational efficiency. 

After a brief review of the IRP method for the solution of the attitude evolution equa- 
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tion, the angular acceleration kinematic model is presented. Applying minimum mean square 
error (MMSE) estimation theory to the perturbation model, the measurement processing al- 
gorithm is developed for both GPS carrier phase signals and vector observations. An attitude 
matrix orthogonalization procedure, incorporated to enhance the algorithm’s accuracy and 
robustness, is then introduced, followed by a derivation of the prediction stage. Two numer- 
ical examples axe then presented, which demonstrate the performance of the new algorithm. 
Concluding remarks are offered in the last section. 

II Integrated-Rate Parameters 

Consider the matrix differential equation 

V(t) = w(t)V(t), HW = V 0 (1) 

where V(t) £ R"", W(t) = - W T (t ) for all t > to, = / and the overdot indicates the 
temporal derivative. Defining 


A{t,t 0 )± fw{T)dT 
JtQ 

Wo(t)±W(t)-(t-t Q )W(t) 


( 2 ) 

(3) 


it can be shown 11 that the following matrix- valued function is a third-order approximation 
of V(t ): 

V{t, to) & {/ + A(t, to) + 4!M) + + t^Js [„(t, to)Wo(t) - Wo(t)A(t, to)] } K„ 

(4) 
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Moreover, V is a third-order approximation of an orthogonal matrix, i.e., V(t, t 0 )V T (t, t 0 ) = 
I + 0 ((£ — t 0 ) 4 ) where 0(x) denotes a function of x that has the property that 0{x)/x is 
bounded as x — ► 0. 

In the 3-D case, the off-diagonal entries of A(t, t 0 ), termed integrated-rate parameters, 
have a simple geometric interpretation: they are the angles resulting from a temporal- 
integration of the three components of the angular velocity vector 

OJ2(t) U*(t)F (5) 

where u>i is the angular velocity component along the i-axis of the initial coordinate system, 
and i = 1,2,3 for x,y,z, respectively. The orthogonal matrix differential equation (1) is 
rewritten, in this case, as 


D(t) = Q(t)D(t), D(to) = D 0 (6) 

where D(t) is the attitude matrix, or the direction cosine matrix (DCM), Q(t) = — M*)xl, 
and [u;(t)x] is the usual cross product matrix corresponding to u>(t). In this case, the matrix 
A(t, to) takes the form 


A(t,to) = -[»(() x] 


where the parameter vector 6(t) is defined as 

^(0 = [0i(O 02(0 03 ( 0 ] 


(7) 


(8) 
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and 


m a f u>i( 

JtQ 


t) dr, z = 1, 2, 3 


(9) 


Let the sampling period be denoted by T = t k+1 - t k . Using the notation 6{k) = 9(t k ), 
the parameter vector at time t k is 9{k) = [^(Jfc) 0 2 (fc) 0 3 (fc)] r and Eq. (9) implies 


9i(k) = [ Ui(r)dT, 
Jto 


i= 1,2,3 


FYom Eq. (10) we have 


9{k 


u/(r) dr 


Define A(A: 4- 1, A) to be the discrete-time analog of A(t, to), i.e., 


( 10 ) 


( 11 ) 


A(k + l,k) = ~[(9(k + 1) - 9(k)) x] 
Also, let ^(k + 1) = - [ip(k + l)x] , where 


( 12 ) 


ip{k + 1) = uj(k + 1) - u(k + 1)T 


(13) 


Then, the corresponding discrete-time equivalent of Eq. (4) is 


D(k + 1) = (/ + A(k + 1 ,k) + l A 2 (k + 1, k + ^A 3 (fc + 1, k) 

L 2 ' 6 

+ ^t[a(& + 1, k)<S(k + 1) - ¥(* + 1)A(A: + 1, A;)] |d(A:) 


(14) 
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which, using Eqs. (12) and (13), can be written as 

D(k + 1) = D[0(k + 1) - $(k),u(k + l),u{k + 1), D(k)] (15) 

III Kinematic Motion Model 

To avoid using the uncertain spacecraft dynamic model, the spacecraft angular acceleration 
is modeled as a zero-mean stochastic process with exponential autocorrelation function. The 
acceleration dynamic model is, therefore, the following first-order Markov process, 

cj(t) = —A u>(t) + v(t) (16) 

For simplicity, a decoupled kinematic model is chosen for the three angular rate components, 
i.e., A = diag{Tf 1 ,r 2 _1 ,T 3 _1 }, where {rj}f =1 are the acceleration decorrelation times associ- 
ated with the corresponding body axes. The driving noise is a zero-mean white process, with 
power spectral density (PSD) matrix 

Q{t) = 2AE 2 , E = diagfaq, a 2 , a 3 } (17) 

The noise variances in Eq. (17) were chosen according to the Singer angular acceleration 
probabilistic model, 9 in which the angular acceleration components, {u;j}f =1 , can be 1) equal 
to u>Mi with probability pM i} 2) equal to —d>\n with probability p^f i , 3) equal to zero with 
probability p 0i , or 4) uniformly distributed over the interval [-%„%)] with the remaining 
probability mass. Using this model, it follows that 

= ^d(l + 4p Mt -p 0i ) (18) 
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The parameters a >*/,, p^, and po, are considered as filter tuning parameters. As customarily 
done, they are selected by experience with real and simulated data, so as to optimally adapt 
the filter to the characteristics of the problem at hand. 

Defining now the system’s state vector as 

*(<) = [t 9 T {t ) u T {t ) <Ju T (t)] T (19) 


the state equation is 



o 

o 


0 

x(t) = Fx(t) 4 v{t) = 

0 0 I 

x(t) + 

0 


0 0 -A 


u(t) 


( 20 ) 


with obvious definitions of F and v(t). Corresponding to the sampling interval T, the 
discrete-time state equation is 


x(k + 1) = <&(T)x(k) 4 v(k ) 


( 21 ) 


where the transition matrix is 


$(T) = = 


I TI A _2 (e _Ar — I + TA) 
0 / A -1 (/ — e -AT ) 


0 0 


,-AT 


and v(k ) is a zero-mean, white noise sequence, with covariance matrix 


( 22 ) 


Q{k) A f e^dtagfO, 0, Q(t) }e f ' r < T -‘> dt (23) 

Jo 


Explicit computation of the integrals in Eq. (23) yields the following expressions for the 
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Fig. 1: GPS Phase Difference Measurement Geometry 
entries of the symmetric covariance matrix Q(k ) 


Qu (k) = A“ 4 E 2 (I + 2A T - 2A 2 T 2 + \ A 3 T 3 - e~ 2AT - 4A Te~ AT ) (24a) 

u 

Qn{k) = A -3 E 2 (/ - 2AT + A 2 T 2 - 2e' AT + e~ 2AT + 2 A Te~ AT ) (24b) 

Qn{k) = A" 2 E 2 (I - e~ 2AT - 2 A Te~ AT ) (24c) 

Q 22 (k) = A -2 E 2 (4e -AT - 3/ - e~ 2AT + 2AT) (24d) 

Q 23 (k) = A _1 E 2 (e _2AT + / - 2e" AT ) (24e) 

Qzs{k) = E 2 (/ — e~ 2AT ) (24f) 


IV Measurement Processing 

GPS Differential Phase Measurements 

Consider the basic GPS antenna array, depicted in Fig. 1. The array consists of the master 
antenna, A m , and the slave antenna, A y These antennas are located on the satellite’s surface, 
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such that the baseline vector between them, resolved in a body-fixed coordinate system, is 
bj. It is assumed that the entire system consists of m b antennas, in addition to the master 
antenna, so that there exist m b independent baselines. It is also assumed that at time £ fc+l , 
m, GPS satellites axe in view. 

Consider the ith satellite, and denote the sightline (unit) direction vector to that satel- 
lite, resolved in an inertial coordinate system, by s*. Let D(k + 1) be the attitude matrix 
transforming vectors in the inertial coordinate system to their body-fixed system represen- 
tations at time tk+i- Let Nij(k + 1) and A <f>ij{k + 1) denote the integer and fractional paxts, 
respectively, of the phase difference between the two carrier signals, corresponding to the ith 
satellite, as acquired by the antennas A m and A,. Denoting by A the GPS carrier wavelength, 
the true (noiseless) signals satisfy 


[A^(* + 1) + N^k + 1)] A = bjD(k + l)sj (25) 

The standard GPS carrier wavelength is 19.03 cm. In th is work, it is assumed that the integer 
part of the phase difference between the two receivers is known from a previous solution. 1,12 

In practice, the phase measurements will be contaminated by noise, the primary source 
of which is due to the multipath effect. 1 Denoting the noise corresponding to the baseline bj 
and the sightline s* by n^k + 1), the real measurement equation is 

[A 4>ij{k + 1) + Nij(k + 1)] A = bjD(k + l)sj + hij(k + 1) (26) 

where it is assumed that n^k + 1) ~ N(0, a 2 (A: + 1)). Typically it can be assumed that the 
noise standard deviation is on the order of 5 mm. 1 From Eq. (26) we obtain the normalized 
measurement equation 


A 4>ij{k + 1) + N tJ (k + l)=bj D(k + l)s, + n tJ {k + 1) (27) 
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where we have defined bj = bj/X and n t j(k+ 1) = h tJ (k + 1)/A. The normalized measurement 
noise satisfies n t] (k + 1) ~ K(0, a^(k + 1)), where 


&ij(k +1) — a tj (k + 1)/A. 


(28) 


GPS Measurement Linearization 

At ijt+i the minimum mean square error (MMSE) predicted vector is x(k + l\k), and its 
corresponding prediction error covariance matrix is P(k + 1| A:) = E{x(k + l|/c)x r (A; + l|fc)}, 
where the estimation error is 


z O' 1*0 - x(j) ~ x(j\k). 


(29) 


Using Eq. (15), Eq. (27) is rewritten as 

Nij(k + 1) + &4>ij{k + 1) = bjD[d(k + 1) - 6(k),u(k + l),ti;(A: + 1), D(A:)]sj + Uij(k + 1) 

(30) 

Next, we linearize the nonlinear measurement equation (30) about the most recent estimate 
at tk- n, i.e., 



9(k + l|fc) 


S9(k+ 1) 

x(k + 1) = x (k + l|fc) + Sx(k + 1) = 

ui{k + 1 1 A:) 

+ 

Sui(k + 1) 


Cj[k "h l\k) 


6uj(k + 1) 


where 69(k + 1), Su>(k + 1) and 5<l>(k + 1) are the perturbations of the state components 
about the nominal (i.e., predicted) state. Let D*(k\k) denote the a posteriori, orthogonahzed 
estimate of the attitude matrix at time t k , to be discussed in the next section. Using now 
the most recent estimates for D(k) and x(k ), namely £>*(A;|A:) and x(k\k), respectively, in 
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Eq. (30), it follows that 


A 0 y (* + 1) + Nij(k + 1) = b]D[6{k + l|fc) + 89(k + 1) - 9{k\k),u{k + 1|A:) + 6u(k + 1), 

u(k + 1|A) + Su>(k 4- 1), + Tiij(k + 1) (32) 


As discussed in the sequel, the a posteriori IRP estimate is zeroed after each measurement 
update (due to full reset control of the IRP state). We will, therefore, use the reset value of 
the IRP estimate, 9 c (k\k) = 0, in Eq. (32). Now expand D about the nominal state using a 
first-order Taylor series expansion, i.e., 


D^9{k + 1 1 At) *+* 89{k +• 1), tI/(A -f- 1|A) -1- 8u>{k -1-1), u>(k 1|A) -h 8ui(k -f- 1), ( Ar | A;) J 

= D{k + nfc) + ± ± mM± ± ii*»w*)] I mk+1) 


»=i 




dD[d(k + l\k)Mk + l),u(k + 1[A), £>*(A;|A:)] 


t=l 

3 




du>i 

dD[9{k + l\k), u{k + l|Jfc),cj(A: + 1), D*(k A:)] 


Ul(fc + l|fc) 


«=i 


du>i 


u(k+ l|fc) 


0(fc4*l|fc) 
6cJi(k 4 “ 1 ) 

+ 1 ) 


(33) 


where (•) denotes ‘evaluated at £’ and 


D(k + 1 1 A:) ^ D[9{k + l\k), u(k + 1|£),<1>(A: + l|Jt), D*(k\k)] (34) 


Differentiating Eq. (14), the sensitivity matrices appearing in Eq. (33) are computed as 

p\ 

QQ-D[9{k + 1), u{k + 1|A), t5(A + l|Ar), D*(k\k)] = Gi[9(k + 1), ip(k + 1|A)]D*(A;|A;) (35a) 

O -I 

D[9(k + l\k),u(k + 1 ),dj{k + 1 1 A:) , D'{k\k)] = -TFi[9(k -1- l|Jfc)] D m (k\k) (35b) 

C/U/j O 

r) 1 

— D[0(fc + l|A:),w(* + 1| k),u(k + 1), D'{k\k)] = -~T 2 F t [9{k + l|A)]D*(ifc|A) (35c) 

OtJi (} 
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for i — 1,2, 3, where 


4>{k + l\k) ± u(k + l\k) - Tut(k + l\k) (36) 

and 

Gi(0,rp) = ^(9ef + e t 0 T ) - 6J - (1 - ^ ||(9|| 2 )[e i x] + - e^ 7- ) + ^0i[0x] (37a) 

Fi{9) = ei e T - 9e[ (37b) 

where ei is the unit vector on the ith axis, i — 1,2, 3 . 

Using Eqs. (33), (35) and (37) in Eq. (32) yields 

A 4>ij{k + 1) + Nij(k + 1) — bjD(k + l|fc)sj = hjj[k + l)5x(k + 1) + u ij(k + 1) (38) 

where the observation vector hij(k + 1) € R 9 is defined as 


hij(k + 1) = ^e^(/c + 1) hJlj(k + 1) + l)j (39) 

and the elements of the vectors heiffk + 1) G R 3 , h ui j(k + 1) e R 3 and h^ffk + 1) e R 3 are 

h, 0 ijp (k + 1 ) = bjG p [9(k + l\k),ip(k + l\k)]D*(k\k)si, p=l,2,3 (40a) 

hu ijp (k+ 1) = ±TbjF p [9(k + l\k)]D*(k\k) Si , p= 1,2,3 (40b) 

hcjij p (k + 1) = — Th ui j p (k + 1), P = 1, 2, 3 (40c) 

Define now the effective GPS measurement to be 

yfjffk + 1) = A 4> l3 {k + 1) + Nij - bj D(k + 1| k) Si (41) 
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Then, using this definition in Eq. (38) yields the follov/ing scalar measurement equation: 


yfj{k + 1) = h xj T (k + 1 )Sx(k + I) + n XJ (k -I- 1) (42) 

For the m b baselines and m 3 sightlines, there exist m s x m b scalar measurements like Eq. (42). 
We next aggregate all of these equations into a single vector equation, such that the mea- 
surement associated with the baseline bj and sightline s, corresponds to the pth component 
of the vector measurement equation, where p — (j — 1 )m 3 + i. This yields 

y <t> (k + 1) = H^ik -I- 1 )Sx(k + 1) + n <t> { k + 1) (43) 

where the pth row of the matrix H <t> (k -I- 1) is hij T (k + 1), the measurement noise satisfies 

n?(k + 1) ~ N(0, R*(k f 1)) (44) 

and the covariance i?^(A: + 1) is a diagonal matrix whose diagonal elements are 

J£(* +!) = »§ (45) 


Vector Observation Aiding 

If the sole source of attitude information is the GPS carrier phase signals, then Eq. (43) 
should serve as the basis for the development of the m<asurement update algorithm (in the 
next section). In the case that vector observations are available, this information structure 
needs to be augmented. 

Assume that a new pair of corresponding noisy vect )r measurements is acquired at tk+i- 
This pair consists of the unit vectors u(k -f 1) and v(k + 1), which represent the values of 
the same vector r(k -f 1), as modeled in the reference coordinate system and measured in 
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the body coordinate system, respectively. The direction-cosine matrix D(k + 1) transforms 
the true vector representation u 0 into its corresponding true representation v 0 according to 

v 0 (k + 1) = D(k + l)uo(* + 1) (46) 


Ass um ing no constraint on the measurement noise direction, the body-frame measured unit 
vector, v(k + 1), is related to the true vector according to 


v(k + 1) = 


Vo (k + 1) + n(,(fc + 1) 

IM&+ 1) +n' v ( k + 1)11 


(47) 


where the white sensor measurement noise is n' v (k + 1) ~ >f(0, R„(k + 1)). Since both 
v 0 (k + 1) and v{k -I- 1) axe unit vectors, it follows from Eq. (47) that 


v(k + 1) = v 0 (k + 1) + n v (k + 1) 


(48) 


where n„(A;+l) = T4(^+ 1 ) n U fc + 1 ) and y i( fc + 1 ) - I~v 0 (k+l)v$(k+ 1) is the orthogonal 
projector onto the orthogonal complement of span{vo(A; + 1)}. To a good approximation, 
the effective measurement noise is a zero mean, white Gaussian sequence with covariance 


Rvik + 1) = ?i(k + l)K(k + l)K( k + !) ( 49 ) 


To account for non-ideal effects (e.g., star catalog errors), it is assumed that the modeled 
reference vector is related to the true vector according to 


u{k + 1) = Uo(k + 1) + n u (k + 1) 


(50) 


where n u ± u 0 is a zero mean, white Gaussian noise, that is uncorrelated with n„ and has a 
known covariance matrix Ru(k). 
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Vector Measurement Linearization 

Using Eq. (15), Eq. (46) can be rewritten as 

v 0 (k+ 1) = D[d{ k+ 1) -6(k),u(k + l),u(k + l),D(A;)]iio(fc + 1) (51) 

Linearizing about the predicted estimates and using Eqs. (31), (48) and (50), it follows that 

v(k + 1) - n v (k + 1) = D[6(k + 1|&) + 86{k + 1 ),u(k + 1|A:) + 8uj(k + 1), 

u>(k + l|fc) + 5u>(k + 1), .D*(A;|A;)] [u(fc + 1) — n u (k + 1)] (52) 

where, as previously done in the GPS measurement linearization, the reset value of the IRP 
estimate, 9 c (k\k) = 0, has been used. Expanding C about the nominal state using the 
first-order Taylor series (33) yields 

3 

v{k + 1) - D(k + l|A:)u(A: + 1) = ^ [(?* [0(A; + 1|A:),^(A: + l|fc)]^(A: + 1) 

1=1 L 

+ ^TFi[6(k + ll*;)]^*; + 1) - ^T 2 F t [§(k + ll*)]^* + 1 )]d*(A:|A:) U (A: + 1) 

— D(k -I- l|A:)n u (A: + 1) 4- n v (k + 1) 

= H v (k + l)5x(k + 1) - D(k + l|fc)n u (A; + 1) + n v (k + 1) (53) 

where the observation matrix H L/ (k + 1) is written in block matrix form as 

H v (k + 1) = [//i(A; + 1) H 2 (k + 1) H 3 (k + l)j (54) 
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and the columns of the submatrices H^k + 1) e R 3 - 3 , i = 1, 2, 3 are 


H\j(k + 1 ) = G } [8(k + 1|A:),^(A: + 1|A:)]£>*(A:|A;MA; + 1) (55a) 

H 2j (k + 1) = ^TFjleik+llk^D'ikltyuik + l) (55b) 

Hjj(k + 1) = —T H 2 j(k + 1) (55c) 

for j = 1, 2, 3. Notice that the same sensitivity matrices are used here, as in the linearized 

GPS measurement equation, which implies obvious computational saving. Define now the 
effective measurement and measurement noise to be, respectively, 

y u {k + 1) ^ v(k + 1) - D(k + l|A:)it(fc + 1) (56) 

n v {k + 1) = n v (k + 1) — D(k + l|fc)n u (A; + 1) (57) 

Then, using these definitions in Eq. (53) yields the following measurement equation: 

y v (k + 1) = H u (k + 1 )5x{k + 1) + n v (k + 1) (58) 

where n v {k + 1) ~ X (0, R w {k + 1)) is the white measurement noise, and 

R v (k + 1) 4 ^(k + 1) + D{ k + l|jfc)J2„(A + 1 )D T {k + l|Jfc) (59) 


Measurement Update 

To process the measurements, define now 


A 

y = 

y* 

. H = 

H* 

A 

, n = 

n* 


y v 


h r 


n v 


( 60 ) 
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where n ~ N(0, R ) and R = diag{i? 0 , R u }. Since 

5x(k + 1) = x(k + 1) - x(k + l|fe) = x(k + 1| k) (61) 

and x(k + l|fc) is an unbiased, MMSE predictor, we have 

E{8x(k + 1)} = E{x(k + 1|A:)} = 0 (62) 

and 

cov {Sx(k + 1)} = cov{x(fc + 1 1 Ac) } = P(k + 1| k) (63) 

thus 

5x(k + 1) ~ X(0, P(k 4 1| k)) (64) 

Using the linearized measurement equation and the stai istical properties of the measurement 
and prediction errors, the MMSE estimator of 5x(k + 1) is 

6x(k + l\k + l) = K(k + l)y(k + l) (65) 

where K(k 4- 1), the estimator gain matrix, is computed as 

K(k + 1) = P(k + l\k)H T (k + 1) [H(k + l)P(k + l\k)H T (k + 1) + R(k + 1)] _1 (66) 
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Also, Sx{k + l|fc + 1) = x{k + 1| A: + 1) - x{k + 1|A) which, used in Eq. (65), yields the state 
measurement update equation 

x(k + l|Jfe + 1) = x(k + l\k) + K(k + 1 )y(k + 1) (67) 

Subtracting x(k + 1) from both sides of the last equation yields 

x(k + 1| k + 1) = [/ - K(k + 1 )H{k 4- 1 )}x(k + 1|A:) - K{k + 1 )n(k + 1) (68) 

from which the resulting covariance update equation is 

P(k + l|Jfc + 1) = [/ - K{k + \)H{k + 1)] P(k + 1|A) [/ - K(k + 1 )H(k + 1)] T 

+ K(k+ l)R(k + l)K T (k + 1) (69) 

where the filtering error covariance is P{k + 1|£ + 1) — E{x(k + 1 1 A: + l)x T (k + 1|A: + 1)}. 

To compute the measurement-updated attitude matrix at time tjt+i, we use the most 
recent estimate x(k + l\k + 1) and the estimated attitude matrix corresponding to time t k 
in Eq. (14). This yields 

E){k + l\k + 1) = |/ + A(k + 1, k) + \^ 2 {k + 1, k) + gi 3 (£ + 1, k) 

+ ^r[A(k + 1, k)$!(k + 1|A: + 1) - 4>{k + l|fc + 1 )A(k + 1, A:)] J D*{k\k) (70) 

where the a posteriori estimates of A(k + 1, k) and ^(k + 1) are defined, respectively, as 

A(k + 1, A) = — [0(& + 1|A + l)x] (71) 

$(Jfc + l|* + l)4-[^(Jfc+l|* + l)x] (72) 
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where 


ij>(k + l|fc + 1) = u(k + 1 1 A: + 1) - Tdj{k + 1|A + 1) (73) 

and D*(k\k) is the a posteriori, orthogonalized estimate of the attitude matrix at time t k , to 
be discussed in the next section. 

Finally, since the a posteriori attitude matrix, D(k + l|fc + 1), is computed based on 
the a posteriori estimate, 9(k + 1|£ + 1), this implies a full reset control [13, p. 332] of the 
parameter vector, i.e., 


e c (k+l) = d(k + l)-6(k~l\k + l) (74) 

where 9 c (k 4- 1) is the reset state vector at t k + 1 , anc. a corresponding reset of the state 
estimate, 

6 c (k+l\k + l) = 0 (75) 

which is then used in the ensuing time propagation step. Since the reset control is applied 
to both the state vector and its estimate, no changes are necessary in the estimation error 
covariance matrix. 


V Attitude Matrix Orthogonalization 

To improve the algorithm’s accuracy and enhance its stability, an additional orthogonaliza- 
tion procedure is introduced into the estimator, following the measurement update stage. In 
this procedure, the attitude matrix closest to the filtered attitude matrix is computed. 
Given the filtered attitude matrix D(k + 1 1 A: + 1), the attitude matrix orthogonalization 
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problem is to find the matrix 


D*(k + l\k + 1) = arg min D(k + 1|A: + 1) — D (76a) 

oeR 3 ’ 3 

subject to 

D T D = I and detZ? = +l (76b) 

Being a special case of the orthogonal Procrustes problem, 14 the matrix orthogonalization 
problem can be easily solved using the singular value decomposition (SVD). 15 Thus, if 

D(k + 1|A + 1) = U(k + 1)E(* + 1) V T (fc + 1) (77) 

is the SVD of the matrix D(k + l|fc + 1) where U(k + 1) and V(k + 1) are the left and right 
singular vector matrices, respectively, and T,(k + 1) = diagjsi, S2, S3} is the singular value 
matrix where Si > S2 > S3, then 

D*{k + l\k + 1) = U(k + l)diag{l, 1, detU(k + 1) det V(fc + l)}V T (k + 1) (78) 

In real-time attitude determination and control the excessive computational burden as- 
sociated with the SVD might render its use prohibitive. In such cases, the following ap- 
proximate orthogonalization method, consisting of a single-step application of the iterative 
method introduced in Ref. 16, can be utilized: 

D*{k + 1|* + 1) = N(k + l)D(k + l|Jfc + 1) (79) 
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where 


N{k + 1) = ^/ — ^D(k + l|fc + l)D T (k + 1|A: + 1) (80) 

Remark 1. Using an approach similar to that used in Ref. 17, it can be shown that, to 
first-order, the orthogonalization procedure does not affect the statistical properties of the 
estimator and, therefore, does not necessitate any adjustments in the algorithm. 

VI Prediction 

In the prediction step at the reset a posteriori estimate at time t*, x c (k\k) (computed 
with the reset IRP estimate) and its corresponding error covariance matrix, P(k\k), are 
propagated to time tk+i- 
Using Eq. (21), we have 

x(A:+l|A:) = $(r)x c (c|fc) (81) 

Using this result with Eq. (21) yields the covariance propagation equation 

P(k + l|fc) = $(T)P(fc|A:)3> T (T) + T(T)Q(k)T T (T) (82) 

To propagate the attitude matrix to 1 we use the most recent IRP, attitude-rate and 
angular acceleration estimates, and the orthogonalized DCM estimate corresponding to 
in Eq. (14). This yields 

D(k+l\k) = \l + A(k + 1 , lb) 4- l-A 2 (k + 1 , k) + \A 3 {k + 1, k ) 

^ 2 6 

+ A{k + 1, k)V{k + 1| k) - V(k + 1| k)A(k + 1, fc)] |D*(ifc|Jfc) (83) 
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where the a priori estimates of A(k + 1, k) and ^(k + 1) are defined, respectively, as 


A(k + l,k) = ~[6(k+ l|jfe)x] 

(84) 

V(k + l\k) = -[^(fc+l|*)x] 

(85) 


VII Numerical Study 

Two numerical examples axe presented in this section, to demonstrate the performance of 
the new estimator, and illustrate the performance enhancement achieved by using aiding 
vector observations. 

Example I 

In this example, three non-orthogonal baselines were used: bi = [l.O, 1.0, 0.0 ] T , & 2 — 
[0.0, 1.0, 0.0] T , 63 = [0.0, 0.0, 1.0 ] T . Two fixed sightlines were observed at all times, 
Si = ^-[l.O, 1.0, 1.0] T and s 2 = ^[0.0, 1.0, 1.0] r . The non-normalized GPS signal noise 
standard deviation was 5.0 mm. When vector measurements were used, the noise equivalent 
angle of the inertially-referenced observations was set to 5.0 arc-s, while the- body-referenced 
vector measurements were simulated to be acquired by a low accuracy attitude sensor with a 
noise equivalent angle of 0. 1 deg. These measurements corresponded to a randomly selected 
vector, which was kept constant throughout the run. 

The angular rates of the satellite satisfied u>i{t ) = AiSin(jrt + fa), where the amplitudes 
Ai are 0.02, 0.05 and 0.03 deg/s, the phases fa are 7r/4, 7 t /2 and 37 t /4 rad, and the periods 
T t are 85, 45 and 65 s for i = 1, 2, 3, respectively. The initial angular rate estimates were 
all set to zero. The true initial attitude corresponded to Euler angles of 30 deg, 20 deg and 
10 deg in roll, pitch and yaw, respectively, while the filter’s initial state corresponded to 
Euler angles of 25 deg, 15 deg and 5 deg, respectively. The filter was run at a rate of 20 Hz, 
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and the measurement processing rate was 10 Hz. The Singer angular acceleration model was 
used with parameters set to r = 10 s, u> M = 10 -4 rad/s 2 , p M = p 0 = .001 for all three axes. 

In Fig. 2, the three true Euler angles are shown for a typical run. These angles were 
computed from the true attitude matrix assuming a 3-2-1 angle sequence. The Euler angle 



Thm <•> Tim* {•) 


(b) Pitch angle 


(c) Yaw angle 


Fig. 2: Example I: True Euler angles. 


estimation errors, computed from the estimated attitude matrix assuming a 3-2-1 Euler angle 
sequence, are shown in Fig. 3. Fig. 4 presents the angular rates estimation errors for the 
same run, with and without vector measurement aiding. The mean and standard deviation 
of the estimation errors are summarized in Table 1, which demonstrates the effect of the 
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on Euler angle estima- 












Table 1: Example I: The effect of vector observation aiding on estimation per- 
formance. 



GPS only 

Vector observation aiding 

Mean 

Standard 

deviation 

Mean 

Standard 

deviation 

Roll angle (deg) 

-9.2 x 10" 4 

5.8 x 10" 2 

-4.2 x 10" 3 

3.9 x 10- 2 

Pitch angle (deg) 

3.0 x 10- 3 

8.1 x 10- 2 

— 1.4 x 10~ 3 

2.2 x 10" 2 

Yaw angle (deg) 

7.1 x 10~ 3 

9.5 x 10" 2 

9.9 x 10~ 4 

2.2 x 10- 2 

tui (deg/s) 

-5.8 x 10~ 4 

9.5 x 10~ 3 

5.1 x 10“ 5 

8.2 x 10- 3 

u/ 2 (deg/s) 

2.0 x 10- 4 

2.7 x 10“ 2 

-2.8 x 10" 4 

1.2 x 10- 2 

oj 3 (deg/s) 

4.0 x 10~ 5 

1.5 x 10~ 2 

3.9 x 10" 4 

6.5 x 10" 3 


vector observation aiding in reducing the estimation errors standard deviation. 

Example II 

In this example, the same parameters were used as in Example I, except for the follow- 
ing. The three baselines used were now = [0.1, 1.0, 0.l] T , 62 = [0.0, 1.0, 0.0] T , 63 = 
[0.0, 0.0, 1.0 ] 1 . As can be observed, the first two baselines are almost colinear. The angular 
velocity of the satellite was u> = [0, 236, 0] T deg/hr, which is typical for an Earth-pointing, 
low Earth orbit satellite, with pitch rate of one revolution per orbit. The Singer angular 
acceleration model parameters were set to r = 10 s, ujm = 10 -5 rad/s 2 , Pm — Po = 001 for 
all three axes. As in the first example, vector measurements, when available, corresponded 
to a randomly selected vector, which was kept constant throughout the run. 

In Fig. 5, the three true Euler angles are shown for a typical run. Fig. 6 shows the Euler 
angle estimation errors (the estimated angles were computed from the estimated attitude 
matrix assuming a 3-2-1 Euler angle sequence). Fig. 7 presents the angular rates estimation 
errors for the same run, with and without vector measurement aiding. The estimation error 
statistics are presented in Table 2. As can be observed, especially from Fig. 6 and Table 2, 
the robustifying effect of aiding the GPS measurements with vector observations is very 
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Table 2: Example II: The effect of vector observation aiding on estimation per- 
formance. 



GPS only 

Vector observation aiding 

Mean 

Standard 

deviation 

Mean 

Standard 

deviation 

Roll angle (deg) 

7.2 x 10" 3 

6.4 x 10~ 2 

1.3 x 1(H 

2.0 x 10~ 2 

Pitch angle (deg) 

1.1 x 1CT 3 

3.8 x 10“ 2 

-4.8 x lO" 4 

2.0 x 10" 2 

Yaw angle (deg) 

7.7 x 10~ 3 

8.7 x 1CT 2 

4.6 x 10~ 3 

2.2 x 10" 2 

uji (deg/s) 

-9.6 x 10" 6 

4.8 x 1CT 4 

2.5 x 10" 6 

3.3 x 10~ 4 

u) 2 (deg/s) 

2.8 x 1CT 5 

9.9 x 10- 4 

3.7 x lO’ 5 

5.1 x 10~ 4 

u 3 (deg/s) 

3.4 x 10" 6 

9.3 x 10" 4 

-5.8 x 10" 6 

3.5 x 10" 4 


significant in this ill-conditioned case. 

VIII Conclusions 

A nonlinear sequential estimator has been presented, that uses differential GPS carrier phase 
measurements to estimate both the attitude matrix and the angular velocity of a spacecraft. 
The algorithm is based on the IRP third-order minimal parametrization of the attitude 
matrix, which is at the heart of its computational efficiency. Avoiding the use of the typically 
uncertain (and frequently unknown) spacecraft dynamic model, the filter uses a polynomial 
state space model, in which the spacecraft angular acceleration is modeled as an exponentially 
autocorrelated stochastic process. When vector observations are available (e.g., from low 
accuracy Sun sensors or magnetometers), the estimator’s structure can be easily modified 
to exploit this additional information and, thereby, significantly enhance the algorithm’s 
robustness and accuracy. Numerical examples have been presented, that demonstrate the 
performance of the proposed algorithm and the advantages of aiding the GPS carrier phase 
signals with vector observations, even when the vector measurements are of relatively low 
accuracy. 
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